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Diffusion processes with branching play an important role in statistical dynamics. They are a 
common approach to the computing of quantum mechanical groundstates, p TJ and serve as models 
for population dynamics and as physical pictures for biological evolution, n|. 

On a computer the efficiency of this simulation method is limited by the approach to the infinites- 
imal time step, which is necessary to perform alternating diffusion and branching steps. 

In this paper, a method is described, which eliminates the infinitesimal time step for a certain class 
of branching processes, if the process of interest can be "embedded" into another process, which is 
solvable by other analytic and/or numerical methods. The simplest choice for the embbeding process 
is given by a process with a constant branching rate, which dominates the rate of the embedded 
process. 
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I. INTRODUCTION 

A population of random-walkers (j)(x,t) which inde- 
pendently move, multiply and die shall be described by 
an equation of motion where the movement is given by 
a diffusion process, while the branching process depends 
on the spatial coordinate x. The probability for either a 
birth or a death process occuring during an infinitesimal 
time-inter vail [t, t + dt] in the spatial interval [x, x + dx] 
is given by \S(x)\(j>(x,t) d d xdt, where S(x) < denotes 
decay processes and S(x) > denotes birth processes. 
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a birth process occurs during [t, t + dt] in 
the spatial interval [x, x + dx] 

_ ( 0, if S(x) < 0, 

~ \ S(x)<j)(x, t) dt d d x, if S(x) > 0, 

a decay process occurs during [t, t + dt] in 
the spatial interval [x,x + dx] 

-S(x)<j)(x,t)dtd d x, i£S(x)<0, 
0, if S(x) > 0. 
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To stress the close relation to imaginary time quantum- 
mechanics the reproduction operator S is introduced as 
the negative potential S(x) = —V(x). The simplest case, 
where the the diffusion is homogeneous and isotropic the 
dynamics is governed by the equation 



",t) = {i€-V(x)}</>(x,t) 



(5) 



This equation is equivalent to an imaginary time Schro- 
DiNGER-equation, with the dynamics 



- H = \dl - V(x) 



Or 



-H(j>. 



(G) 



The formal analogy between imaginary time Schrodin- 
GER-equation and statistical physics has been studied 



since |l| , a contemporary approach is given in @. Most of 
the works focuses on the relation between the imaginary 
time Schrodinger— equation and equilibrium statistics. 
However, this paper emphasises, that there is a more nat- 
ural interpretation as a non-equilibrium process, which 
has been explicitely stated in Q and which is implicitely 
used in i,!. 



II. THEORETICAL FOUNDATION FOR THE 
IMAGINARY TIME Schrodinger-EQUATION 



In this section, the method to measure the energy of 
the quantum-mechanical groundstate and to stabilize the 
population shall be derived. In contrast to the common 
bilinear product (<j)\H\(j)), a flux through the population 
shall be used to measure the energy. The population is 
stabilized on the average, by subtracting the groundstate 
energy Eq, which is determined "on the run". 



A. Formal solution 



In analogy to real time Schrodinger— equation 
— = H(f> equation (|^) can be solved using the same 
eigen -function system of H with Hip n = E n ip. This 
way one obtains exponentially growing (decaying) con- 
tributions for Ei < (Ei > 0) , instead of oscillating ones 



(x,t) 



tpn{x), 



(<hWn). (7) 
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Of course, the groundstate ip n grows faster (or decays 
slower) than all other states and will dominate the ensem- 
ble for t — > oo, thus on the long run the normalised popu- 
lation converges to the quantum-mechanical groundstate 



tl>o(x) — lim </>(x,t) 



(x',t)dx' 



(8) 



B. Measuring the energy 

It is remarkable, that the energy of tj> is related to a 
flux through the population, which is defined on an area 
Gby 



F(t) = / H(f>dx 

>G 



<)G 



Vcj)d d x. (9) 



Obviously, an exploding or decaying system cannot be 
in equilibrium. Of course, even if the total population 
is stable J G <f>(x, t) ^constant, there are areas, where the 
random-walkers are born, and other areas reached by dif- 
fusion, where they die. Thus, the RDS cannot be related 
to an equilibrium system. However, under certain condi- 
tions it is convenient to have a stationary flux, which is 
equivalent to a system with a stable population. 



C. Stabilizing the population 

The main idea to stabilize the population is to de- 
termine the energy of the groundstate and to subtract 
it from the effective dynamics. This has the effect of 
putting the lowest eigen-value to zero by hand. As 4> is 
the formal solution of 



0(t) =exp (-#*)<£(()), 



(10) 



it is useful to notice, that for an arbitrary time- 
depending function E{t) the equation 

d t ct>(x,t) = {E(t) - H} <t>(x,t), (11) 

with the same c n as defined in equation (Q), is solved by 

cl>{x,t)=Y,Cneti E{t ' )dt '- Eit i> n {x). (12) 

n 

Thus, 4> can also be stabilised by subtracting an average 
energy with the property J Q E(t') dt' = Eot, for t — > oo, 
which is automatically implied if lirnt_,oo E(t) = E Q . 
Therefore, the average 



E(t) = / F(t') dt' 



JG 



(x',t')dx' dt' 



(13) 



is chosen converging E(t) — ► Eq because <f> — > ifjo f° r 
t — > oo by fll) and inserting <fi into definition (^) . 



III. SIMULATION WITHOUT INFINITESIMAL 
TIME-STEPS 

The starting point is the common time-step algorithm. 
The simulation needs N time-steps to evolve from time 
t = to the final time T, each step of length e = T/N. 
With Rnd and £ denoting standard uniform and normal 
random-numbers the simplest formulation is: 

for n <- to N 

1. t <— TIE, 

2. move x r <— x r + sfe^, 

3. if \S(x r )\e > Rnd then 

perform a birth (S > 0) or death (S < 0) process, 

4. next n. 

At first, it is assumed, that there is an upper bound 
S > S(x) allowing to split step || into two parts: 

3a. if \S\e < Rnd reject any branching, next n, 
3b. if \S(x r )\ < Rnd S then 

perform a birth (S > 0) or death (S < 0) process. 

The probability of a branching event is the probabil- 
ity p a — \S\e to pass step 3a times the probability 

Pb = \S(x r )\ J S to pass step 3b without rejection. The 

result is the correct branching rate p a Pb = \S(x r )\e. 

Now there is an inner loop, l-3a, and an outer loop, 
l-3b. In the inner loop there is no need to evaluate the 
exact branching potential S(x). The probability, not to 
leave the inner loop after one step is Q 1 = (1 — Se). The 
probability, not to leave the inner loop after n steps, or 
equivalent until time t = ne is given by Q n = (1 — Se) n — 
q e (t). The limit e — > can be evaluated 

q e{t ) = {l-Se)^ s -^{e-' s ) t = e-' st . (14) 

Thus the whole loop can be replaced by: 
while t < T 

1. draw an exponentially distributed random-number 
t <- log(l - Rnd) /S, t<-t+T 

2. move x r <— x r + y/r£, 

3. if \S(x r )\ < Rnd S then 

perform a birth (S > 0) or death (S < 0) process, 

4. next exponential time step. 

This way, the time r between to proposed branchings 
is related to the upper bound S for the potential |5(a;)|, 
if the branching has been succeeded or rejected. The rate 
for evaluating the true S (x) is S, which can be chosen to 
be of a typical order of magnitude for the given problem. 
According to the considerations to stabilise the ground- 
state, section (II C), we choose S(x,t) = E(t) — V(x). 

Although, the proof of the algorithm has not been 
given for time dependent S(x,t), this generalisation is 
obvious. Furthermore, in the long time limit, the fluctu- 
ations of E(t) fade away. 
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IV. THE REPRESENTATION OF THE WAVE 
FUNCTION 

The simulations have been performed by ensembles of 
R idependent random- walkers, each described by two 
variables, the coordinate £^ and the time 6„, where 
and when the r-th walker has been a "candidate" for 
its n-th branching process, independent of a branch- 
ing has been performed or not. Since the last branch- 
ing proposal at time at location the random-walker 
has done a free random-walk. When it is considered 
to branch again at time t, its new position is given 
by a normal distributed random-number with mean 
£n and variance a = Wt — Q T n - With g(x,a 2 ) — 

(27rer 2 ) d ^ 2 exp (— x 2 /2a 2 )the wave-function is repre- 
sented by 

R 

<Kx,t)=^g(x-e r ,t-6r) (15) 

r=l 

and averaging the flux (^) reduces to the computation 
and summation of integrals 

f F(t')dt' = j2Yl w-^-i.C-i) 

Jo r=1 „=i 

0^<t 

T 

with U(t, $) = J dt' J d d x V(x) g (x - £, t') . (16) 
o 



V. EXAMPLES 

The quality of the simulations depends on the compu- 
tation of U (t, x). For all standard potentials like polyno- 
mials and COULOMB-Iike potentials, including interact- 
ing electronic systems, U can be solved analytically. 

A. Polynomials and the harmonic oscillator 

In Jl6| ) the function U (r, x) has been defined as the 
time average of the energy. For an arbitrary polynomial 

p ( x ) = T,k u ...,k d 9ki,...,k d Ilti x f onc finds 

d ki 

u(t, X )= ?* *„IIX>*.* J (it) 

fei,...,fed i=lji=0 

with Kjjs = 1+( 2 " 1)J (*) {^2) ■ { j/2+i which is eas y to 
implement. For the d-dimensional harmonic oscillator 
V(x) = P(x) = ^\x\ 2 , one obtains the simple result 

U(t,x) = it\x Q \ 2 + ^t 2 . (18) 



(x r }±5x r 




time t 

FIG. 1. paths of the harmonic oscillator for a single ran- 
dom-walker, x r ± s/t — &n to demonstrate the variance of the 
normal distribution. The diffusion is free for a time r of the 
order of a typical time scale. 

B. The CouLOMB-potential 

Although the Coulomb \j\x\ singularity has no up- 
per bound, the function U can also be computated by 
analytical methods. 

By scaling the HAMlLTONian gets H = — %d 2 + a/\x\ 
with the SOMMERFELD-constant a being the only free 
parameter, thus Eq = -i« 2 . The integrals ( |l6| ) are 
solved in polar coordinates performing the integrations 
in the order ip — > 8 — > r — » t, with the angles defined 
according to the axis — > x. The /Tin-function reads 

U H (t,x) =aV^i u(\x\ A/2i) , (19) 

with u(x') = (l/(2x')+x') erf{x') + e~ x ' 2 /tt - x'. The 
function u(x') has a very simple structure, and is tabu- 
lated at the beginning of each run. The integrals evalu- 
ated for the computation of the CouLOMB-potential may 
easily be generalised to simulate the ortho-helium atom. 
The HAMlLTONian 

2 

H =-J2(^ d l - 2c */M) + a/\xi - x 2 \ (20) 

i=l 

allows a separation of the 6 + 1-dimensional integration. 
For Vi — 2a/\xi\ there is one integration for Xi and an- 
other one for the coordinate which Vi is independent of. 
The V12 — a/\xi — x^\ integration may be transformed 
into the centre of mass system, resulting in a constant fac- 
tor and another integration for the relative coordinates 
which reduces to a CouLOMB-integration again. The fi- 
nal result reads 

2 

U He (t, (x u x 2 )) = 2J2U H (t, Xi) + U H (t, x x - an). (21) 

i=l 
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VI. THE QUALITY OF THE RESULTS 

To keep the number of random-walkers constant, ev- 
ery time a random-walker multiplies (dies) another one 
dies (multiplies) with probability l/(R— 1). For this "an- 
tagonist" the actualization of the coordinates has to be 
carried out, too. This way, R does no fluctuate any more, 
otherwise the fluctuations of R would cause a diffusion 
process of R, which is avoided for convenience. 
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FIG. 2. Systematic error of Eg for the hydrogen atom, de- 
pending on the number of random-walkers. The straight line 
represents 1/7?. 

The antagonist rule introduces an unknown additional 
coupling among two walkers considered as independent. 
This results in a law, the large number of R— 1 pos- 
sible partners weakens the individual coupling. Thus the 
order of the systematic error is higher than the statistical 
error of i?~ 1//2 and can be neglected, if R is sufficiently 
large (Fig. ||). Typical values are 0.1 percentage error 
for R = 2 10 = 1024. 

For harmonic or Coulomb potentials using a constant 
S(x,t) the condition S > S(x) is violated. However, due 
to the fast decay of the wave function this error decays 
exponentially fast for the harmonic oscillator. For the 
Coulomb like problems it only affects a neglectible area. 
Fig. H gives an impression of two typical runs for the he- 
lium problem approaching the experimental groundstate 
energy better than 0.1%. 



VII. CONCLUSION 

The algorithm has proved its capability to simulate 
imaginary time Schrodinger— equations very efficient 
and accurate. The method of embedding has been ex- 
tended to all problems, which can be formulated as 
branching processes, and it has been generalised to ar- 
bitrary time depending rates S(x, t) M. 
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FIG. 3. Eo for typical run of a helium simulation, comput- 
ing time vs. percentage error 

Without respecting the bosonic or fermionic properties 
of the wave function, in multi particle simulations there 
is a mixture of symmetric and anti-symmetric compo- 
nents. Because the bosonic wave-function has a lower 
energy than the fermionic one, for long times the solu- 
tion converges to the bosonic groundstate. Thus, the 
algorithm in its present implementation is ideal for the 
simulation of interacting bosons |7|] . 

One interesting property of the algorithm for 
quantum-mechanics is the fact, that Coulomb 2- 
particle integrals can be solved analytically. Its impor- 
tance is limited by the question, if it will be possible to ex- 
tend it to fermionic systems. For all kinds of Coulomb- 
systems an identity like (^l)) holds as well and the singu- 
larities may be integrated and smoothed out. At least, 
it should be possible to introduce the nodes of a related 
problem as absorbing boundaries similar to 

However, the main advantage is opposite to ||Q that 
there is no need of a "guiding wave function" . Therefore, 
the algorithm is easier to implement. The related ideas 
concerning fermionic systems and excited states [||, 
could be applied as well. 
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